function power = for_max(in)
%优化计算不同分布情况下的功率power
%in代表输入的参数

mir_wide=in(1);%镜宽
mir_length=in(2);%镜长
mir_hight=in(3);%镜高

%% 已知参数
month=9;
D=daysact(datetime(2023,3,21),datetime(2023,month,21));
D=D+365*(D<0);%日期
t=9:1.5:15;%时刻
lat=39.4/180*pi; %纬度
H=3; %海拔 km
delta=asin(sin(2*pi*D/365)*sin(2*pi*23.45/360));%赤纬角
tower=[0,0,80];%吸收塔坐标

%% 镜子离散化
dx=20;dy=20;
[xx,yy]=meshgrid(-mir_wide/2:mir_wide/dx:mir_wide/2,-mir_length/2:mir_length/dy:mir_length/2);
[X,Y]=distribution(mir_wide);%计算各个镜子的坐标位置
mir_center=[X,Y,mir_hight*ones(size(X))];%计算镜子中心坐标

dHR=tower-mir_center;
V_out=normr(dHR);%光的反射向量
dHR=sqrt(sum(dHR.^2,2));%计算镜子到吸收塔的直线距离

%% 计算
for temp_t=3
    ST=t(temp_t);
    omega=pi/12*(ST-12);
    alpha_s=asin(cos(delta)*cos(lat)*cos(omega)+sin(delta)*sin(lat)); % 太阳高度角
    gamma_s=real(acos((sin(delta)-sin(alpha_s)*sin(lat))/cos(alpha_s)/cos(lat))); % 太阳方位角
    gamma_s=(2*pi-gamma_s)*(ST>12)+gamma_s*(ST<=12);
    
    V_in=[cos(alpha_s)*sin(gamma_s),cos(alpha_s)*cos(gamma_s),sin(alpha_s)];%光的入射向量
    V_mid=normr(V_in+V_out);%通过入射光向量与反射光向量计算镜子法向向量
    angle_X=acos(V_mid(:,3));
    angle_Z=angle(V_mid(:,1)+1i*V_mid(:,2))+pi/2;
    
    %各个效率计算
    eff_sb=zeros(size(X));
    eff_cos=V_mid*V_in';
    eff_at=0.99321-0.0001176*dHR+1.97e-8*dHR.^2;
    eff_trun=0.92;
    eff_ref=0.92;
    
    dni=1.366*(0.4237-0.00821*(6-H)^2+(0.5055+0.00595*(6.5-H)^2)*exp(-(0.2711+0.01858*(2.5-H)^2)/sin(alpha_s)));
    
    for i=1:length(X)
        H1_i=[xx(:),yy(:),zeros(length(xx(:)),1)];
        matrix_change_i=martix_trans(angle_X(i),angle_Z(i));
        H2_i=(matrix_change_i*H1_i')'+mir_center(i,:);
        difference=([X,Y]-[X(i),Y(i)]);
        distan=abs(difference(:,1)+1i*difference(:,2));%记录镜子之间的相互距离
        distan_sort=sort(distan); % 从小到大排序
        close_6=find(distan<=distan_sort(6) & distan>0);%选择距离最近的6个镜子
        shadow_index=[];%用于记录阴影遮挡镜子序号
        %% 寻找遮挡镜子
        for j=1:length(close_6) 
            norm_j=V_mid(close_6(j),:); % 镜的法矢量
            jiao_s=(dot(norm_j,mir_center(close_6(j),:))-H2_i*norm_j')/dot(norm_j,V_in); % 这是i镜上各个离散点沿着光射入的方向去找与j镜所在平面的交点参数s
            
            node_loc=H2_i+jiao_s*V_in; % 交点的全局坐标    
            trans_matrix_j=martix_trans(angle_X(close_6(j)),angle_Z(close_6(j))); %  % 给出j镜的旋转变换矩阵，
            node_in_j=(trans_matrix_j\(node_loc-mir_center(close_6(j),:))')'; % 求逆,给出交点在j镜中的局部坐标
            shadow_index=[shadow_index;setdiff(find(abs(node_in_j(:,1))<mir_wide/2 & abs(node_in_j(:,2))<mir_length/2 & jiao_s>0),shadow_index)]; % 判断交点是否在j内,不重复记录已经判定为遮挡交点 
        end
        num_shadow=length(shadow_index); % i镜上的离散点总的被遮挡个数,计算遮挡效率
        eff_sb(i)=1-num_shadow/(dx+1)/(dy+1);%更新阴影遮挡效率
    end
    
    %% 计算光学效率与热输出功率
    eff=eff_at.*eff_cos.*eff_ref.*eff_sb.*eff_trun;%光学效率
    power_out=dni*mir_wide*mir_length*sum(eff);%热功率
    day_power=power_out;
end
    power=-mean(day_power);%日均热功率 由于求最大值 所以需要加负号
end

